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ABSTRACT 

In two previous papers dSalvador-Sole et al.lf2012allbh . it was shown that: i) the typical 
structural and kinematic properties of haloes in (bottom-up) hierarchical cosmologies 
endowed with random Gaussian density perturbations of dissipationless collisionless 
dark matter emerge naturally from the typical properties of peaks in the primordial 
density field and ii) halo statistics are well described by the peak formalism. In the 
present paper, we use these results to model halo substructure. Specifically, making 
use of the peak formalism and the fact that accreting haloes evolve from the inside-out, 
we derive the subhalo mass abundance and number density profile per infinitesimal 
mass for subhaloes of different masses, as a function of the subhalo maximum circular 
velocity or mass, before and after the tidal truncation of subhaloes by the host po- 
tential well. The subhalo properties obtained by assuming that subhaloes are mainly 
made of diffuse particles are in very good agreement with those found in current high- 
resolution iV-body simulations. We also predict the subhalo properties in the opposite 
extreme case, likely better suited for the real universe in CDM cosmologies, that haloes 
are made of subhaloes within subhaloes at all scales. 
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1 INTRODUCTION 

Halo substructure is attracting much attention for its impli- 
cation s in th e problem of galaxy formation. iKlypin et al.l 
1 19991 ) and iMoore et al l \ 19991) noticed that numeri- 
cal simulations of cold dark matter (CDM) cosmolo- 
gies found many more subhaloes at the galactic scale 
than observed. Although the situation has notably 
changed si nce this finding from the observational view- 
point (e.g. IWilliams fc McKed Il997l; iBelokurov et al.ll2006l : 



iKoposov et alJl201Cl : iBullock et al.ll2010h . the uncertainty 
remains on the consistency of CDM cosmol ogies with the 
observed abundance of dwarf galaxies (e.g. IStrigari et al.l 
l20ld : lLovell et all 120121 ). Moreover, the distribution of sub- 
halo maximum circular velocities found in simulations has 
also been recently shown to be in conflict with that ob- 
served in the Milky Way satellites (|Bovlan-Kolchin et al.l 
l201ll : IVera-Ciro et al.|[201ll ). 

The accurate determination of subhalo properties has 
been hampered for a long time by the extremely high dy- 
namic range required by A^-body simulations addressing this 
issue. In the last decade, however, there has b een an impres- 
sive i mprovement (e.g. Ichigna et al.l 19981 : Springel et al 



20011: iHelmi et al.l 12002: Stoehr et alTl2003l: iDiemand et al 



2004bl : lKravtsov et a l. 2004; De Lucia et al.ll2004l ; lGao et al 
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120041: iReed et al.l|2005l). The sim ulations bv lDiemand et all 
(|2007l . I200SD and lSpringel et all (|2008bl . hereafter SWV) fi- 
nally converged to a well-determined subhalo mass function 
for Milky Way mass systems in the concordance cosmol- 
ogy. More recent simulations with similar resolutions have 
begun to study subhalo properties in host h aloes with dif- 
feren t masses, redshifts and concentrations (jAngulo et al.l 
120091 : lElahi et al.ll2009l : IKlypin et al1l201ll ; ICao et al.ll201ll ). 

The situation is less satisfactory, however, from the 
theoretical viewpoint. The origin of substructure is well- 
understood: it is the consequence of the halo hierarchical 
growth and the fact that, when haloes are captured by a 
more massive object, they are not fully destroyed by the 
host tidal field. But the details of the process are not clear 
enough and the models developed s o far (e.g. iFuiita et al.l 
l2002i : IShetbll2003l : lOguri fc Leell2004l ) are still unable to re- 
cover the results of A-body simulations. This is embarrass- 
ing not only because of the need to fully understand this 
phenomenon, but also because, due to their high cost in 
CPU time, current A'-body simulations with still rather lim- 
ited dynamic range inform us on substructure in haloes with 
only a few masses and redshifts. 

In two recent papers, ISalvador-Sole et alj (|2012al : 
l2012bl . hereafter SVMS and SSMG, respectively) have shown 
that the smooth structural and kinematic properties of 
haloes (related to the 1-particle probability distribution 
function) endowed with dissipationless collisionless dark 
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matter emerge from the properties of peaks in the density 
field at an arbitrarily small cosmic time, determined by the 
power-spectrum of density perturbations in the (bottom-up) 
hierarchical cosmology under consideration. This allowed 
these authors to justify the peak formalism ( Doroshkevichl 
ll97d ; lDoroshkevich fc Shandarinll 19781; iBardeen et alj|l986t 
hereafter BBKS; iPeacock fe Heavens! 1 199(A iBond fe Myers! 



Il99rj). more specifically, the ri g orous versio n of it developed 
by iManrique fc Salvador-Solel ill995l; 1 19961 hereafter MSSa 
and MSSb; see also lManrique et al.lll998h and show that this 
is a very useful tool to deal with halo statistics. 

In the present paper, we use this formalism and the 
model for the smooth structural and kinematic properties of 
haloes developed in SVMS and SSMG to describe the prop- 
erties of CDM halo substructure (related to the n-particle 
probability distribution function, with n ^ 2). Although 
iV-body simulations are nowadays able to study clumps-in- 
clumps (e.g. SWV), the only complete studies carried out 
so far on the properties of halo substructure concern first- 
level subhaloes. For this reason, we focus on this particu- 
lar kind of substructure, although we also draw conclusions 
on their internal structure. Specifically, we derive the sub- 
halo abundance and number density profile per infinitesimal 
mass, as a function of subhalo maximum circular velocity 
and mass, before and after the tidal truncation of subhaloes 
by the host potential well. Although the formalism is devel- 
oped under the assumption that haloes form by pure accre- 
tion (PA), the results obtained are shown to be also valid 
for haloes having undergone major mergers. The theoretical 
predictions obtained for the CDM concordance model with 
(fi m , Oa, h, as) — (0.3,0.7,0.7,0.9) are in very good agree- 
ment with the results of numerical simulations. 

The paper is organised as follows. In Section [2j we re- 
cap the main results obtained in SVMS, SSMG, MSSa and 
MSSb in connection with the present study. These results 
are used in Section [3] and Section [4] to derive the sub- 
halo abundance as a function of maximum circular veloc- 
ity and the number density of subhaloes per infinitesimal 
original (non-truncated) mass, respectively. The correction 
for tidal truncation is given in Section \S\ In Section [5] we 
discuss the effects of dynamical friction. Our results are 
summarised in Section [7] A package with the numerical 
codes used in the present paper is publicly available from 
www . am . ub . es/ ~cosmo/haloesftpeaks . tgz. 



2 THEORETICAL BASIS 

As shown in SVMS, in accretion periods between consecutive 
major mergers haloes develop outwardly by keeping their 
instantaneous inner region unaltered. As a consequence of 
this growth, the radius r encompassing the mass M in a 
(triaxial) virialised halo grown by PA (i.e. having suffered no 
major merger) from a (triaxial) protoobject with outward- 
decreasing spherically averaged density profile and Hubble 
flow-dominated kinematics at an arbitrarily small time U 
(where the protohalo is in linear regime) satisfies the relation 



3 r\£ p (M)\ 
10 GM 2 



= 1 . 



(1) 



where £ P (M) is the total energy within the sphere with mass 
M of the spherically averaged protohalo. Equation |l} allows 



one to 
hence, 
In 

0(r) = 



infer the mass profile M(r) of the virialised halo and, 
its spherically averaged density profile (p)(r). 
SSMG, it was shown that the relation 



= 1 



T t V) 



(2) 



(3) 



between the anisotropy profile, /3(r), and the rms scaled po- 
tential fluctuation profile, ((<5$/(<I>)) 2 ) 1/ ' 2 (r), related to the 
same profile in the protohalo, ((<5$ p / ($ p )) 2 ) 1/ ^ 2 (r p ), and the 
generalised Jeans equation for anisotropic triaxial systems 
(see SVMS), 




d[(p)a 2 /(3-2/9)] , 2/3(r) <p)(r 



dr 



+ - 



r 3 - 2/3(r) 



(J) 



for the velocity variance, cr 2 (r), allow one to derive the /3(r) 
and <r(r) profiles for haloes from the triaxial shape of peaks 
acting as their putative seeds by PA. 

The typical properties of peaks, namely the quantities 
£ P (M) and ((<5$ p / ("I>p}) 2 ) 1/2 (rp), in the primordial den- 
sity field were derived from the power-spectrum of density 
perturbations making use of the peak formalism. This for- 
malism is based on the peak Ansatz stating that there is 
a one-to-one correspondence between haloes and peaks as 
suggested by PA. Specifically, any halo with M at the time 
t is traced by a peak in the density field at U filtered by 
means of a Gaussian window, with density contrast <5 p k and 
filtering scale R{, given by 



tfpk(t) = <fc(t) 



cm 

G(t) 



R f {M) = - 



3M 



1/3 



(5) 



where G(t) is the cosmic growth factor, q is the radius, in 
units of _Rf, of the collapsing cloud with volume equal to 
M over the mean cosmic density at ti, p(ti), and 5 c (t) is 
the critical linearly extrapolated density contrast for cur- 
rent collapse. For the concordance model, the best values of 
gand 8 c (z) are respectively 2.75 and 1.93+ (5.92 - 0.472z + 
0.0546z 2 )/(l + 0.000568z 3 ), where z is the redshift corre- 
sponding to the cosmic time t. 

As haloes accrete, their associated peaks describe con- 
tinuous trajectories in the <5 p k--Rf diagram. (The peaks over 
the trajectory are connected between each other in the sense 
explained in SSMa; see also SVMS.) The typical peak tra- 
jectory leading to a halo with Mo at to is the solution of the 
differential equation, 



= —x c (6 P k, Rt) a2(Rf)Rf , 



(6) 



for the boundary condition 5 p k(£o) at Rf(Mo), where 
£c(<5 P k, -Rf) is the inverse of the average (close to the most 
probable) inverse curvature x (equal to minus the Laplacian 
over 02), 



(Ri, 5 p k) = 



(277) 



-1/2 



(1 - 7 2)l/2 



da; - f(x) 



(7) 
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in peaks with <5 p k and R{ (BBKS), being 



/(*) = ^ < ovf 



2 

UttJ 



sx 1 / 2 



(I) 



+ erf 



UJ 2 



31x 2 8\ 
— + 5 6 



5 6 



(8) 



where 7 and a;* are respectively defined as af/(aoa2) and 
75 p k/o"o in terms of the j-th order spectral moments for the 
power-spectrum P(k), 



<T?(Rt) = 



(9) 



The first equality in equation © relates the derivative of 
S p k(t) to the typical mass accretion rate, r^ BB (M,t), of 
haloes with M at t (see MSSb) . 

When a halo suffers a major merger, the 5 p k(iif) trajec- 
tory of its associated peak is interrupted (there is no peak at 
the immediately larger scale to be connected with). At the 
same time, a new peak appears (there is no peak at the con- 
tiguous smaller scale to be connected with), with the same 
density contrast as the disappeared peak but at a substan- 
tially larger scale, that traces the new halo formed in the 
major merger. Major mergers are the only way peak trajec- 
tories are interrupted. When a halo is accreted by another 
much more massive halo, the trajectory of the associated 
peak is not interrupted, so the peak becomes nested into 
the collapsing cloud of the larger scale peak with identical 
density contrast tracing the accreting halo. On the other 
hand, when a peak disappears because the associated halo 
suffers a major merger, its trajectory is interrupted but not 
those of its nested peaks, which also become nested in the 
collapsing cloud of the new peak resulting from the merger. 
In this way, a complex system of peak nesting at multiple 
levels is built similar to the nesting of subhaloes in haloes 
(see Sec. [3j . 

One important implication of the SVMS and SSMG 
model is that the properties of haloes having suffered major 
mergers are indistinguishable from those of haloes grown by 
PA. This result allows one to understand why the typical 
density and kinematic profiles for haloes grown by PA are 
representative of all haloes, regardless of their aggregation 
history, and why the peak Ansatz suggested by PA also holds 
for all haloes, regardless of their actual aggregation history. 

Below, we detail the form of several unconditioned 
and conditional number densities of peaks, nested or non- 
nested within o t her p eaks calculated in MSSa, MSSb and 
iManrique et al.l l|l998l ) that will be used in the following 
sections. Note that all these quantities depend on the power- 
spectrum of the hierarchical cosmology considered through 
the spectral moments defined above. 

The number density of peaks with density contrast <5 p k 
at scales Rf to Rf + dRf is the number density of peaks 
at scale R{ with density contrast 8 greater than 5 p k that 
cross such a density contrast when the scale is increased to 
-Rf + dRf or, equivalently, with 8 satisfying the condition 



<5 P k < 8 < <5 pk + cr 2 (-Rf) x R t AR { . 



(10) 



Thus, such a number density can be obtained by integrating 
over 8 and x the density of peaks with height v = S/ao(R{) 
and curvature x in infinitesimal ranges, J\f p k(v, x, R{ ) Av Ax, 



calculated by BBKS 

iVpk(-Rf,5pk) d -Rf 

_ (x)(i?f,<5p k ) _!i 0-2 (Rt 



(2ix) 2 Rl 



da; / AvAf p k(v,x,Rf) 



R f AR t 



(T (Rf 



(11) 



where v' pk — ^ p k + [(T2(Rt) / cro(Ri)]RfARf and (a;) is the av- 
erage curvature of peaks with 5 p k and Rf 

<*>(JMpO= ^T^mj dxxf(x)e-^^T . (12) 

Likewise, the conditional number density of peaks with 
8 at scales Rf to Rf + dRf subject to being located in the 
collapsing cloud of a non-nested background peak with 8' at 
R/, N p £ Bt (Ri, 8\Rf', 8') dRf, can be obtained by integrating 
the conditional number density of peaks with 8 at scales Rf 
to Rf + AR{ subject to being located at a distance r, in units 
of qRi, from the background peak, A r p k(Rf , 8\Rf' , 8', r), out 
to q = 1, 

N£ st (Ri,6\R f ',5')=C- 1 [ drSr 2 N ph (Ri,5\Rf',S',r), (13) 
Jo 

with the latter conditional number density of peaks ob- 
tained, as the ordinary number density above, by integrat- 
ing over v and x the conditional density of peaks with 
those arguments in infinitesimal ranges, subject to being 
located at the distance r from a background peak with v' at 
Rf', Npk(v,x,Rf\Rf',6',r)dvdx, calculated by BBKS, with 
8 satisfying the condition l|10j l. 

r°° rKk 

NpkiRfjSpklR/ ,8' ,r)AR{= Ax duM p ^(v, x, Rt\Rc ,6,r), 

JO JupU 

(s)(iM P k,r) - ^"'^'"^ ^(Rf) ( . 

where (x)(Ri, 5 p k> f) is the average curvature of peaks with 
5 p k and Rf at a distance r from a peak, given by 

(27r)- 1/2 f°° -[*-**(*■)]" 



(x)(R[, 5 pk , r) = ■ 



dxa;/(a:)e m-^M] , .(15) 



In equation (|15[) . we have used the following nota- 
tion: e(r) — v/l — e(r) 2 , aj*(r) = j(r)u(r), j 2 (r) — 



l + e(r) 2 { 



2 (i-nr 



n = (Rt/Rhf, Rl = (Rf 2 + i? f ' 2 )/2, 



6(r)=(i? f Rf'/ J R 2 ) (n+3)/2 ;? (r,Rf'), u'(r) = g(r,R f '), 
g(r,R { ') = {l - [A«5'(r)] 2 /a (i?f')} 1/2 and 
~f \ T 1 - r i 



j(r) 1 — e(r) 2 



1 - e(r) 2 ri , , 
''I 1 _ ri ) - e{r)v (r) 



with the mean and rms density contrast at r from a peak, 
8'(r) and [A8'(r)] 2 , respectively equal to 



5'(r) 



■7 V7 



jS ( tp V ip 



a:ao(Rf 
1 -7 2 



+ 50 (Jr.) 



[A8'(r)] 2 = a 2 \l- 



I-7 2 



3*0' V 2 V> 



1p 2 +[ 27^ + — rl 



3(V 



7«^ 



(17) 
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being tj) the ratio £0")/£(0), ip' its r-derivative, £(r) 
the mass correlation function at Rf and u equal to 
(qRt) 2 <J2(Rf) / ao(R{) ■ In equation (|13|) . the factor 

c ' tS;^ > i' dr3r ' w '' w - aia ' i>) ' <i8) 

with s equal to the mean separation between non-nested 
peaks drawn from their mean density (eq. [19] below), is 
to correct for the overcounting of background peaks, as in 
the preceding calculation of N pk (Rf, 5\Rf' , 5', r) they are not 
necessarily non-nested. 

Finally, given the preceding number densities, it is read- 
ily seen that the number density of non-nested peaks with 8 
at scales Rf to R{ + dRf, N(Rf, S)dRf, is the solution of the 
Volterra integral equation correcting the ordinary number 
density of peaks (eq. [TT]) for nesting, 

N(Rt,6) = N pk (Ri,5) 
1 



dR f 'M(R f ') N(R U S) N^ st (R t , S\R/ , S), (19) 



and that the number density of peaks with 5 pk at scales Rf 
to -Rf + dRt , nested in non-nested peaks with S' p]l at scales 
Ri to Rf' + dRf' , is given by 

N ncst (R f -> Rf', S pk -> S' pk )dR { dR f ' 
= N^(Rf, S pk \R { \ S' pk ) dR f jV(R f ', S' pk ) ^2^1 dj R f ', ( 20 ) 

where M(Rf')/p(ti) is the volume of the collapsing cloud of 
the peak with 8' pk at Rf' and N pk st (R t >, S pk \R t ■', S' pk ) is the 
conditional number density of peaks with S pk at Rf subject 
to being located in the collapsing cloud of non-nested peaks 
with 8' pk at Rf' . 



3 SUBHALO ABUNDANCE 

In dark matter clustering, first-level subhaloes develop in 
two ways: i) through the accretion by a halo of much less 
massive partners (with substantially higher concentrations) , 
which become first-level subhaloes of the accreting halo at 
the same time that their own first-level clumps become 
second-level ones and so on; and ii) through major merg- 
ers of similarly massive haloes (with similar concentrations) , 
where the merging objects meld and their respective first- 
level subhaloes are transferred as such to the new halo re- 
sulting from the merger. 

As explained in Section [21 the processes of halo accre- 
tion and major mergers are correctly traced by peak trajec- 
tories in the 5 pk -Rf diagram. Furthermore, the halo-nesting 
they produce is also correctly traced by the corresponding 
peak-nesting. Indeed, when haloes are accreted, the peaks 
tracing them survive and become nested into the collaps- 
ing clouds of those peaks tracing the accreting haloes, while 
the peaks already nested within them become second-level 
nested peaks and so on. On the other hand, in major merg- 
ers, peaks tracing the merging haloes disappear and their 
first-level (second-level,...) nested peaks automatically be- 
come so in the collapsing cloud of the new peak tracing 
the halo formed in the merger. Both behaviours reproduce 
that above mentioned of haloes and subhaloes in accretion 
and major mergers. Thus, by counting the first-level peaks 



nested in the collapsing cloud of peaks in the density field 
at ti, we can estimate the number of first-level subhaloes in 
the associated haloes at t. 

The total number of first-level peaks with 5 pk and scales 
greater than Rs nested within the collapsing cloud of a non- 
nested peak with S pk at scale Rf, N(> Rs, 5 pk \Rf , S pk ), fol- 
lows from equation (|20[l . the result being 

rRf 



N(>Rs,S pk \R f ,5 pk )= 



M 

Wi, 



ARf < Npi (ilf , <Spk | -Rf , <5, 



Rs 

dRf' N^ tncst (Rf,S pk \R { ',S pk ) 



R f 



v P k 



'(Rf', S pk \Rf,S 



pk I 



M(Rf') 

m 



(21) 



In equation (|21[) . the integral over Rf' is to correct 
the number density of peaks nested in the seed of the 
halo for those peaks nested in intermediate-scale peaks 
so as to ensure that only first-level nested peaks are 
counted. The factor N pk st {Rf',S pk \ Rt,S pk ) inside that 
integral comes from the probability for the intermedi- 
ate peaks to be nested in the seed of the halo, equal 
to N™ st (Rf',S pk \Rf,S pk )N(Rf,S pk )/N pk (Rf',S pk ) times 



M/p(ti), and the function N$ ncst (R f , <5 pk |R f ', <5 pk ) 
solution of the Volterra integral equation 



N, 



P k 



(Rf,5 P k|Rf ,S pk ) = N pk s (Rf,<5 P k|Rf ,<5 P k) 
rR{ ' 

dR/' N pk st (Rf, <5 pk |Rf", <5 p k) 

IR, 



7v t 1st nest / 7-j / 

xiV pk (R f 



A IP' A S M W) 

, p k -Kf ,OpkJ Z77T — 

pit;) 



(22) 



gives the conditional number density of peaks with <5 pk at 
scale Rf subject to reach for the first time the same den- 
sity contrast at an intermediate scale Rf' when the scale 
is increased from Rf. This ensures that the correction for 
intermediate nesting is not overcounted. 

As first-level subhaloes with masses greater than Ms in 
a halo with M at t, N(> Ms, t) are correctly traced by first- 
level peaks with <5 pk (i) at scales greater than Rs = Rf(Ms) 
(eqs. [5]) nested in the collapsing cloud of a non- nested 
peak with S pk (t) and Rf(M), their cumulative abundance 
N(> Ms,t) must also be well-estimated by the correspond- 
ing abundance of nested peaks, N(> Rs, 5 pk \Rf, S pk ), given 
by equation (|21|l . In Figure [T] this theoretical cumulative 
subhalo abundance for current Milky Way mass (1.4 x 10 12 
Mq) haloes in the concordance cosmolog y is compared 
to tha t found in numerical simulations by Diemand et al.l 
(|2008h and SWvfl As can be seen, there is excellent agree- 
ment, particularly in the case of SWV results. Figure[5]shows 
the theoretical subhalo mass abundance as a function of the 
scaled maximum circular velocity for different halo masses. 
Except for a small shift at the large mass end, the predicted 
subhalo abundance is essentially independent of halo mass, 
in agreement with a very common idea, though with rather 
limited empirical support. 



1 The halo mass is Mq = 1.4 X 10 12 in all c ases, although 
the mass is defined within r^o in lDiemand et all J2008T) ■ V200 in 
SWV and ~ rgo in the present paper (see above). Nonetheless, 
all three curves overlap at the scale of Figure [T] 
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Figure 1. Theoretical cumulative abundance of non-truncated 
subhaloes (solid red line) for a Milky way mass halo in the concor- 
dance cosmology as a function of subhalo maximum circular ve- 
locity, V m ax, scaled to that of the halo. For compa rison, empirical 
curves obtained by SWV (dashed black line) and lDiemand et al.l 
J2008I 1 (dotted black line). 



As can be seen, the theoretical cumulative subhalo 
abundance shows a small bump at large masses which arises 
from a similar bump in the conditional peak number den- 
sity, (Rf, 6 p k\R{ , 5pk), at scales Rf comparable to Rf. 
This latter function is approximate for Rf close to Rf (see 
iManrique et al,lll99of ): it should vanish when Rf approaches 
Rf more rapidly than it actually doefl This suggests that 
these bumps may be an artefact due to the less steep fall 
of iVpk (Rf, 5 p k|-Rf , <5 p .k) a t large scales. But this is hard to 
ascertain. Empirical data are too noisy there to asses the re- 
ality or not of the bump in the subhalo abundance. In fact, 
there are indications that it is real: had we only slightly 
sanded the bump in the conditional peak number density, 
the resulting subhalo abundance would take negative values. 
For this reason, we have preferred to conserve it and adopt a 
sharp cutoff at Rf /Rf equal to one tenth for subhalo masses 
greater than one hundredth of the host mass and at Rf/Rf 
equal to one hundredth otherwise. Such a cutoff does not 
essentially alter the theoretical subhalo abundance shown in 
Figures [T] and [2] while it notably improves the behaviour of 
the subhalo number density profile derived below for subhalo 
masses close to the host mass. 



2 Not only can there be no peaks nested in other peaks with 
identical scale but also within peaks with slightly larger scale. 
The capture by a halo of another similarly massive one necessarily 
causes a major merger, so the two peaks disappear. 



Figure 2. Same as Figureffjbut for several halo masses (coloured 
lines). 



4 SUBHALO NUMBER DENSITY PROFILE 

Given the inside-out growth of haloes formed by PA (see 
Sec. [2}, the cumulative abundance of subhaloes with masses 
greater than Ms within the radius r, N(< r, > Ms), must 
coincide with the cumulative subhalo abundance by the time 
the halo radius was equal to r. Consequently, the differential 
subhalo abundance, both per infinitesimal halo radius and 
subhalo mass, in a halo with Mo at to is 



N(r,Ms) = \ 



cUV* .(Ms)' 



d<5, 



P k 



d<5, 



5{Rf) 



dRf 



dRf 
dM 



dM 
dr ' 



(23) 



'Rt(M) 



'M(r) 



where 
NsJMs) 



D(M S 



JVpT 4 [Rt (M s ),5 pk \ Rt ( W , V] 



M[ R f (S pk )] j vlM .. ! 
(M s ) 

dRf ' iVpf ncst [R t (Ms), S pk \R { ', <5 pk ] 



Rf(M s ) 



< N pk st [Rf', <5 pk \Rf (<5 pk ), <5 pk ] - 



(24) 



is the differential subhalo abundance obtained by dif- 
ferentiation of the cumulative subhalo abundance 7V(> 
Ms), given in equation (|21l) for M and t respec- 
tively equal to M[i? f (<5 pk )] and t(S pk ), with D(M S ) = 

f 2~l 1//3 

qln [6p 2 (ti)M s ] \ , M(r) given by equation JTJ, R C (M) 

given by equation ([SJ and Rf(S pk ) equal to the inverse typi- 
cal peak trajectory solution of equation ([6]) leading to a halo 
with Mo at to- 
la Figure [3J we show the theoretical number den- 
sity profile per infinitesimal mass for subhaloes with Ms, 
n(r,M s ) = (A-nr 2 )~ 1 N(r, Ms), scaled to the total mean 
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Figure 3. Theoretical non-truncated subhalo number density 
profiles per infinitesimal mass for subhaloes with Mg, scaled to the 
total number, for subhaloes with different masses (solid coloured 
lines), compared to the halo mass density profile of the NFW pro- 
file (solid black line) . The effects of including a primordial diffuse 
particle component that is progressively accreted by haloes is also 
shown (dashed coloured lines). For comparison, we also plot the 
Einasto law fitting the truncated subhalo number density profiles 
in simulated haloes (dotted black line). 



number density of such subhaloes, (n(Ms)}, so obtained. 
This scaled number density profile shows a cutoff, preceded 
by a short bending, at small enough radii that depends on 
the subhalo mass (see e.g. the curves for subhaloes with 
10 8 M Q and 10 6 M @ j3 The reason for this behaviour, also 
found in numerical simulations (see lAngulo et al.l |2009|) , is 
well-understood: there can be no clump with mass Ms inside 
the radius r encompassing that mass because the accretion 
of such a subhalo at the time where the halo had the mass 
M(r) ~ Ms would automatically cause a major merger and 
the consequent destruction of the merging objects. Apart 
from that short bending and cutoff, all the scaled number 
density profiles for subhaloes with different masses overlap 
with the mass density profile of the halo. The ratio between 
the subhalo number density and the total halo mass density 
predicted by the model (see eq. |23|). 



n(r, M s ) 
p(r) 



dS 



P k 




is flat, indeed. This result is at odds with that found in 
simulations, where the scaled number density profiles for 
subhaloes of different masses also overlap with each other, 



but show a much less steep profile than the halo mass density 
profile (see the dotted curve in Fig. [3|. 

The situation in simulations indicates that there is a dif- 
fuse dark matter component outside subhaloes that becomes 
dominant at small radii (SWV). But is this reasonable? In 
idealised hierarchical cosmologies, all the dark matter is ex- 
pected to be locked into virialised haloes of different masses 
that develop through minor and major mergers. And the 
same is true for the matter inside haloes: all of it is expected 
to be locked into subhaloes of different masses. Even if sub- 
haloes are tidally truncated by the host potential well (see 
Sec. [5]), the liberated matter will be in the form of subhaloes 
previously seen as subsubhaloes. It is true that when dark 
matter begins to cluster, after decoupling (or after equality 
if decoupling took place earlier), it is in the form of a dif- 
fuse component which is accreted by the first haloes formed 
by monolithic collapse. But accretion of diffuse matter pro- 
ceeds in a very short time compared to the cosmic times we 
are interested in, so we can neglect such a transient phas^f). 
In simulations, it is instead normal to find some amount of 
diffuse dark matter in current haloes for two reasons: i) sim- 
ulations start with unclustered dark matter at much smaller 
redshifts (of about 100) and ii) haloes (subhaloes) below the 
resolution mass contribute to a melt diffuse component until 
it is fully accreted by more massive haloes (subhaloes). As 
a consequence, only about 60 % of the total mass in cur- 
rent haloes is aggregated through minor and major mergers 
(about 40 and 20 %, respectively); all the rem aining mass 
is acc reted in the form of diffuse dark matter (|Wang et al.l 
l201ll) . This modifies the hierarchical way CDM haloes clus- 
ter at the small mass end and lead s to the presence of di f- 
fuse particles until quite a large z (|Angulo fc Whitc]l2010l ). 
Therefore, if we are to compare the predictions of the model 
with the results of A-body simulations, we must account for 
this effect. 

In the presence of diffuse dark matter, the result above 
that the scaled subhalo number density profile overlaps with 
the halo mass density makes no sense. It would imply that 
the fraction of mass in the form of the diffuse component 
has the same density profile or, equivalently, that the mass 
fraction accreted by haloes in the form of diffuse matter is 
constant, whereas the diffuse dark matter should be progres- 
sively accreted, so its fraction falling into haloes should di- 
minish with increasing time and, given the inside-out growth 
of accreting haloes, with increasing radius in any individual 
halo, causing a downward bending of the subhalo number 
density towards the halo centre as observed in simulated 
haloes. 

To calculate the expected bending of the subhalo num- 
ber density at small radii in a given simulation, we need to 
know the time-evolving mass fraction in the diffuse compo- 
nent outside haloes, /dc(t)- This mass fraction satisfies the 
differential equation 



d/dc 
dt 



fdcjt) 1 
1 - /dc(t) pit) 



dMr^ c {M,t)N{M,t) 



3 The cutoff for subhaloes with 10 10 Mq is located at log(r/R) ~ 
0.1, but it is not preceded by any short bending likely due to the 
effects mentioned above concerning the approximate conditional 
number density of nested peaks at large Rf. 



4 This may not be the case for warm dark matter, whose decou- 
pling marking the beginning of the clustering takes place much 
later. 
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/dc(t) 
l-/dc(*) 



(26) 



where r acc (M, i) is the mass accretion rate of haloes with 
mass M at t, given by equation JBJ, and M a im is the mass 
resolution of the simulation. The solution of equation (|26[) 
for the initial condition /dc(tsim) = 1, with t s im the starting 
time of the simulation, is given by the implicit equation 



/dc(t) - HUc(t)\ = 1 



dtr dc (t) . 



(27) 



Thus, the mass accreted by the halo at t in the form of 
subhaloes is diminished by a factor 1 — fdc{t) compared to 
that in the case of no primordial diffuse component. Given 
the halo inside-out growth, this implies that the contribution 
from subhaloes to the halo mass density profile changes from 
p(r) to p(r){l — fdc[t(r)]}, with t(r) the time where the 
accreting halo reaches the radius r. The effect of such a time 
varying mass fraction in the diffuse component for the initial 
cosmic time i s im corresponding to z = 127 and the resolution 
mass Msim equal to 10 4 Mq as in SWV simulations is shown 
in Figure [3] The curve so obtained is much like the one 
found by SWV, although not identical. As we will see in 
SectionJSJ the difference is likely due to the effects of subhalo 
truncation not considered yet. 

The presence of a diffuse dark matter component should 
have very little effect, however, on the cumulative subhalo 
abundance, N(> Ms), shown in Figures [T] and [2] The rea- 
son is that, despite the outward-decreasing subhalo number 
density profiles, the number of subhaloes increases with ra- 
dius, meaning that they are mostly aggregated by the halo at 
late times when essentially all the diffuse dark matter com- 
ponent has already disappeared (even in simulations) . Only 
if we were analysing the subhalo abundance at very high 
redshifts (or very small radii) should the effect of the diffuse 
dark matter component also be taken into account when 
dealing with the subhalo abundance. Note that the same is 
true for the halo mass function: at very high redshifts it will 
be affected by the diffuse dark matter component, which 
should not be present in the real CDM universe. This is not 
taken into account in studies of that quantity from iV-body 
simulations. 

To sum up, in the case of (essentially) no primordial dif- 
fuse dark matter component, as in the real CDM universe, 
the scaled number density profile for (non-truncated) sub- 
haloes with any given mass should coincide with the scaled 
halo mass density profile. On the contrary, in the case of 
a substantial amount of diffuse dark matter, as in numeri- 
cal simulations or in the real universe soon after decoupling 
(or after the time of equality), the scaled number density 
profiles for subhaloes of different masses should also overlap 
with each other, but not with the halo mass density profile. 
They should be equal to this latter profile times the factor 
1 — /dc(r) giving the mass fraction clustered in haloes by the 
time t(r) when the halo reached the radius r. 

An important consequence of the previous result is that 
the spatial distribution of (non-truncated) subhaloes is the 
same in haloes grown by PA as in haloes having suffered 
major mergers. If it were different in both kinds of haloes, 
then the typical mass density profile would also be differ- 
ent, which would be in contradiction with the results of N- 
body simulations (see SVMS). Strictly, the possibility re- 



mains that the deviation in the typical mass density profile 
for subhaloes of some mass is exactly balanced by that for 
subhaloes of the remaining masses, but such an arrangement 
is very unnatural. We therefore conclude that the spatial 
distribution of subhaloes must be independent of the host 
aggregation history. As discussed in SVMS, this conclusion, 
far from being unexpected, reflects the fact that virialisa- 
tion is a real relaxation process. As such, it must cause the 
memory loss of the initial conditions, not only regarding the 
smooth halo structure and kinematics, but also regarding 
halo substructure (but see Sec. [©J . 



5 THE EFFECTS OF TRUNCATION 

When subhaloes are aggregated by a halo, they are tidally 
truncated by its potential well. Consequently, to compare 
the predictions of the model with the results of numerical 
simulations we must account for the effects of truncation. 
In fact, tidal truncation alters not only the mass of sub- 
haloes but also the number of subhaloes with a given origi- 
nal non-truncated mass, Ms, due to the appearance of new 
first-level subhaloes of that mass, previously seen as subsub- 
haloes, that are liberated from their host subhaloes. 

Let us first concentrate in the change produced in the 
number of subhaloes with a given non-truncated mass. The 
number density per infinitesimal mass of subhaloes with 
non-truncated mass Ms corrected for the number density 
effect owing to truncation is 

N tv (r, M s ) = N(r, Ms) + / dM s N(r, Ms) 

J M min 



R(M S ) 



R tT (r,M s ) 



drN~ (r,M s ). 

M s ,M(r) V ' ° J 



(28) 



where N(r, Ms) is the number density of non-truncated sub- 
haloes with mass Ms, calculated in Section[4] R (r,Ms) is 
the truncation radius of the original subhaloes with mass Ms 
located at r and M™ m is the minimum subhalo mass that 
gives rise by truncation to new first-level subhaloes with Ms. 
The subindexes in the second-level (differential) truncated 

subhalo abundance, N~ (f , Ms), indicate that this sub- 
Ms, A/(r) V aJ ' 

halo number density profile corrected for truncation refers 
to a host, in this case a subhalo, with mass Ms at the time 
when it was aggregated by the halo, with a mass at that 
moment equal to M(r), and hence, different from the mass 
Mo of the halo at to- 

We will consider two extreme cases. In case (a), all 
CDM particles are in subhaloes of a certain mass, as theo- 
retically expected in the real CDM universe at late times, 
so the truncation of first-level subhaloes yields only new 
subhaloes previously seen as subsubhaloes; we thus have 
M™ m = Ms. In case (b), subhaloes are instead essentially 
made of diffuse particles, so the truncation of first-level 
subhaloes does not modify the number of these subhaloes 
(just their mass as well as the total mass of diffuse parti- 
cles in the intrahalo medium); we thus have M™ ln = M. 
Clearly, in case (b), equation (|28[) nas the trivial solution 
A tr (r, Ms) = N(r, Ms), while in case (a) equation ||28) is an 
integral equation neither of Fredholm nor of Volterra type, 
but can still be solved in the way explained in Appendix IA1 



© 0000 RAS, MNRAS 000, 000-000 



8 Salvador-Sole, Serra, & Manrique 




Figure 4. Same as Figure [JJ but for truncated subhaloes in cases 
(a) (dashed red line) and (b) (solid red line) corresponding to sub- 
haloes made of subsubhaloes (and so on) and of a diffuse particle 
component, respectively. The cumulative abundance of truncated 
subhaloes in case (a) coincides with the cumulative abundance of 
non-truncated subhaloes shown in Figure [T] 



The effect of truncation in the cumulative subhalo abun- 
dance as a function of V ma x is shown in Figure H] Note that 
the quantity Vmax is insensitive to the strength of tidal trun- 
cation because the maximum circular velocity in a subhalo is 
reached at a radius smaller than the truncation radius. This 
is the reason why numerical studies usually plot the subhalo 
abundance as a function of Vmax instead of as a function of 
the truncated mass much harder to estimate (see the dis- 
cussion below). In case (b), truncation does not produce 
any apparent change in the cumulative subhalo abundance 
relative to that plotted in Figure [T] recovering that found 
in simulations. The reason is that subhaloes harbour only 
diffuse particles, so the subhalo number does not change 
where subhaloes are truncated. In contrast, we do expect 
an important change in case (a), better suited for the real 
CDM universe, owing to the appearance of new subhaloes 
previously seen as subsubhaloes. As shown in Figure [4] the 
abundance of first-level subhaloes then increases dramati- 
cally (about two orders of magnitude) in comparison with 
the abundance shown in Figure [1] and found in numerical 
simulations. 

All the previous results favouring case (b) indicate not 
only that, in simulations, a large fraction of the mass in 
simulated haloes is in the form of a diffuse component (see 
e.g. SWV), but also that such a diffuse component must 
be widely dominant in subhaloes so that very few new sub- 
haloes emerge by tidal truncation of other subhaloes. This 
does not necessarily mean that there is no subhalo at any 
level higher than one. It just indicates that subsubhaloes 
must be rare enough for not having significant effects on 
the general properties of substructure as drawn from cur- 



Figure 5. Same as Figure|4]but as a function of truncated mass 
Mg r . For comparison, cumulative abundance of non-truncated 
subhaloes (solid black line) as in Figure[T]but as a function of non- 
truncated mass Ms instead of maximum circular velocity V ma x- 

rent high-resolution numerical simulations. SWV report the 
detection of subhaloes up to third level. However, accord- 
ing to the present results, these third-level subhaloes should 
be seen only within the most massive subhaloes and their 
most massive subsubhaloes, so that the total number of sub- 
haloes liberated by truncation would be insignificant com- 
pared to the number of them directly aggregated. There are 
several reasons for such an important lack of subsubhaloes in 
simulations. Subsubhaloes are previously truncated by the 
subhalo potential well and this is also true for third-level 
subhaloes within subsubhaloes themselves and so on. The 
higher the level of subhaloes, the earlier they typically form 
and the less massive they typically are. The earlier subhaloes 
form, the larger is their mass fraction below the mass resolu- 
tion. And, the less massive the subhaloes, the more centrally 
concentrated they are, so the more severe the tidal disrup- 
tion they yield in their own subhaloes. Therefore, we expect 
the mass in simulated subhaloes to be, indeed, mostly in the 
form of diffuse particles. 

But this is not what we expect to find in the real CDM 
universe with (essentially) no primordial diffuse component. 
Neglecting the minimum halo mass, all haloes have grown 
through mergers between less massive progenitors previously 
formed, so substructure should essentially obey case (a). The 
difference between the first-level subhalo abundance pre- 
dicted by the model in cases (a) and (b) implies that there 
should be, in the real CDM universe, two orders of magni- 
tude more first-level subhaloes than usually considered from 
the results of numerical simulations. This, together with the 
fact that in the real CDM universe the subhalo number den- 
sity profile is steeper than found in simulations (see below), 
might have important implications on the detectability of 
CDM from the enhanced flux of cosmic rays produced in its 
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annihilation in n earby subhaloes (e.g. ISpringel et alj|2008bl ; 
lElahi et al.ll2009f ). However, the abundance of dwarf galaxies 
in, say, a Milky Way mass halo is not affected because the 
tidal truncation of subhaloes does not liberate new galaxies 
that were previously hidden. Moreover, luminous (and cold 
baryonic) matter usually lies at the centre of subhaloes, so 
the subhalo abundance relevant for the expected number of 
satellite galaxies rather corresponds to case (b). 

We can now turn to the second effect: the change in the 
mass of the truncated subhaloes. To express the preceding 
subhalo abundance and number density profiles as a func- 
tion of the subhalo truncated mass Mg r [j we must take into 
account the relationship between that mass and the original 
non-truncated mass Ms, 

»fl tr (r,M s ) 



M| r (r,M S ) =4tt 



drr pM s ,M(r)(r), 



(29) 



where PM s ,M(r)( r ) is the typical subhalo density profile 
equal to that for haloes with mass Ms at the time of the sub- 
halo aggregation when the halo had a mass equal to M(r). 

The cumulative abundances of truncated subhaloes for 
Milky Way mass haloes as a function of Mg r are plotted in 
Figure [5] The log-log slopes we find are equal to —1.12 and 
— 1.05 for cases (a) and (b), respectively (or —2.12 and —2.05 
for the differential subhalo abundance). The slope found 
by SWV in their numerical simulations was —0.90 ± 0.03 
(—1.90 ± 0.03), hence once again closer to the value pre- 
dicted by the model in case (b). Note that, although the 
difference between the theoretical and empirical slopes in 
case (b) is small, it may be essential for having a convergent 
or divergent number of subhaloes for masses approaching to 
zero. It is true that this limit is actually not reached due to 
the cutoff in the power-spectrum (and the non-negligible ve- 
locity dispersion) of dark matter particles, but those slopes 
still tell at which extent the mass fraction in low-mass sub- 
haloes is dominant or not. The possible origin of the slight 
departure in the slope between the predictions of the model 
and the results of numerical simulations is discussed below. 
In any case, even if the total number of subhaloes in case (b) 
converged, that in case (a) should diverge as found here, so 
our results point, in the case of CDM cosmologies, to a halo 
mass fraction in the form of low-mass subhaloes (below the 
resolution of current simulations) much larger than usually 
thought. 

The subhalo number density per infinitesimal mass cor- 
rected for truncation as a function of the truncated subhalo 
mass, Mg r , or simply the real truncated subhalo number 
density per infinitesimal mass, takes the form 

/•Al 

N tI (r,M t s r )=N[r,M s (r,M t s I )]+ / dM s N(r,M s ) 

J M™ 



R(M S ) 



R tT (r,M s ) 



M s ,M(r) y ' s ; 



(30) 



with the function Ms(r, Mg r ) implicitly defined by equa- 

5 In simulations, the subhalo mass is usually taken equal to the 
truncated mass Mg r plus the unbound mass belonging to the halo 
in the volume occupied by the subhalo, denoted by M su ]->. We have 
checked that the use of M su ]-> instead of Mg r does not significantly 
alter the results presented here. 



tion (|29[) . The integral equation (|30[) can be solved for 
N tT (r, Mg r ) in the two extreme cases (a) and (b) above in 
the same way as equation 1)281) for N r (r, Ms). Note that the 
truncation radius R tT in equations (|28|) and (|30|) depends, 
for a given halo mass, not only on the radius of the subhalo 
at the aggregation time and its non-truncated mass, but also 
on its orbit, which in turn depends on the host kinematics, 
modelled in SSMG. 

The truncation radius is very hard to determine in nu- 
merical simulations. This is why different authors adopt dif- 
ferent pr ocedures leading to so mewhat different results. For 
instance. iDiemand et al. (|2007l ) took the radius at which the 
density of the subhalo (corrected for the local background 
contribution) is equal to the local background density, while 
SWV adopted the radius at which the clump mean inner 
density (also corrected for the local background contribu- 
tion) is 0.02 times the mean inner background density, shown 
to lead to a truncation radius in ove rall agreement with 
the th eoretical tidal radius defined by Bin nev &: Tremaind 
l|1987l ) . In the present mode l, we adopt the better m o tivate d 
truncation radius given in iGonzalez-Casado et al.l l|l994f ). 
These authors showed that, regardless of the shape of the 
orbit, subhaloes are truncated by the host potential well es- 
sentially at the radius encompassing an inner mean density 
equal to that of the host halo at the clump pericentre. Note 
that such a truncation radius would roughly coincide with 
Diemand et al.'s provided clumps described circular orbits; 
unfortunately, this is not the case in general. On the other 
hand, it would coincide with SWV. truncation radius pro- 
vided the halo mean inner density at the subhalo pericentre 
were 0.02 times the local halo density, which is, in general, 
not the case either. Assuming subhaloes with the median 
velocity for a normal distribution with radial and tangen- 
tial velocity dispersions given by the model in SSMG, we 
determined the typical pericentre reached by subhaloes lo- 
cated at any given radius. Then, assuming the non-truncated 
subhaloes at their aggregation time with the typical halo 
density profile with the mass - conce ntration (M-c) relation 
given by ISalvador-Sole et all l|2007l . hereafter SMGE0), we 
calculated, from the halo mean inner density at the resulting 
subhalo pericentre, the wanted subhalo truncation radius. 

Figure [6] shows the theoretical number density pro- 
file per infinitesimal mass for truncated subhaloes with 
M s r , scaled to the total number of such subhaloes for 
Milky Way mass haloes in cases (a) and (b) compared to 
those found in numerical simulations or, more exactly, to 
the shallow profile of the Einasto form fitting them. In 
numerical simulations, these scaled number density pro- 
files are indeed fo und to be much shallower than the 
halo density profile |l5iemand et al. 200 4bl;lGao et al.ll2004l : 
iNagai fc Kravtsovll2005l ; lDiemand et al.ll2007l; SWV) and in- 
dependent of subh alo mass i Diemand et al.l l2004al ; SWV; 
iLudlow et ai1 l20091. In the left panel of Figure[6l we see that 
the truncated subhalo number density profiles predicted by 
the model in case (a), with no primordial diffuse compo- 
nent, are on the contrary similarly steep as the NFW profile 
fitting the profiles found for non-truncated subhaloes and 
show now a slight dependence on subhalo mass. The sit- 



6 The M-c relation provided by SMGH is consistent with the 
SVMS model for CDM haloes (see SVMS). 
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Figure 6. Same as Figure[3]but for truncated subhalo masses Mg r , in the same cases (a) (dashed lines) and (b) (solid lines) as in previous 
figures, compared now to the Einasto law (solid black line) that fits the subhalo number density profiles drawn from simulations (SWV). 
For comparison with Figure [3] we also plot the NFW law fitting the predicted number density of non-truncated subhaloes (dotted black 
line). Left panel: predictions for the case of no primordial diffuse particle component. Right panel: case (b) predictions for the case of a 
primordial diffuse particle component that is progressively accreted by haloes. 
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Figure 7. Same as right panel of Figure l6l but using the lZhao et al. ( 2009) (left panel) and iKlypin et all J201ll) (right panel) M-c toy 
models instead of the SMGH physical model to calculate the effects of subhalo truncation by the host potential well. 
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uation does not improve in case (b): the number density 
profiles show a more marked dependence on subhalo mass 
and are much shallower than found in simulations at large 
radii, while they become steeper at small radii. And an in- 
termediate case between (a) and (b) would not improve the 
results: the theoretical number density profiles would always 
show a clear dependence on A/g r and be convex instead of 
concave. 

As mentioned, the fact that the empirical subhalo num- 
ber density profiles are independent of subhalo mass and 
shallower than the halo mass density profile implies that, in 
numerical simulations, a substantial fraction of dark matter 
is in the form of a diffuse component that increases inwards. 
According to the discussion in Section [4] a bending of the 
theoretical number density profiles towards those found in 
simulations is expected, indeed, in the case that there is a 
primordial diffuse component that is progressively accreted 
by haloes. Such an effect was calculated in Section U for 
non-truncated haloes. Thus, we have repeated the same cal- 
culations for truncated haloes in case (b) (in case (a) there 
is no diffuse component). That is, we have considered that 
the contribution of subhaloes to the halo mass density is 
given by p(r)[l — /dc(f)] instead of p(r). Then, the log-log 
slope of the truncated subhalo abundance does not essen- 
tially change, but the truncated subhalo number density 
profiles do markedly. As shown in the right panel of Fig- 
ure |S] they then become essentially in agreement with the 
results of numerical simulations. 

The small deviations that still remain (in the curva- 
ture of the profiles and in the dependence on subhalo mass) 
are likely due to the SMGH M-c relation used in the mod- 
elling of truncation for very large z, which apparently is 
not accurate enough. To see the kind of effect the adop- 
tion of one particular M-c relation has on these results, 
we have repeated the same calcu lation with two diff erent 
M-c relations: th ose provided by iKlypin et al.l l|201lf) and 
IZhao et alj (|20091 ). The predicted log-log slopes of the differ- 
ential truncated subhalo abundance are then equal to —2.07 
and —2.01 using Klypin et al. M-c relation and —2.04 and 
—2.01 using Zhao et al. M-c relation for cases (a) and (b), 
respectively (to be compared with the slopes of —2.12 and 
—2.05 found using the SMGH M-c relation and the slope 
of —1.9 found in numerical simulations by SWV). The trun- 
cated subhah2juimber_density profiles in case (b ) predicted 
using IKlypin et all l|201ll ) and lZhao etafl (|2009T ) M-c rela- 
tions are shown in Figure [7] As can be seen, the theoretical 
subhalo number density profiles so obtained deviate more 
markedly from the profiles drawn from simulations and are 
more mass-dependent than those shown in the right panel of 
Figure [6] Therefore, the predictions drawn from the SMGH 
M-c relation are neatly preferable. It may seem strange that 
the M-c relations drawn from numerical simulations give 
poorer results than th ose derived fr o m th e SM GH model. 
We no te, however, that IKlypin et all <|201ll ) and lZhao et al.l 
(2009) M-c relations do not actually fit the results of simu- 
lations; they are the extrapolations of the real fitting expres- 
sions, through some guessed toy models, to the much wider 
mass and z domains involved in the present calculations. In 
any event, these results clearly show that any slight devia- 
tion of the M-c from the true relation has notable effects, 
indeed, in the theoretical truncated subhalo number density 
profile. 



Before concluding this Section, it is important to re- 
mark that, as the non-truncated subhalo number density 
profile, N(r, Mg 1 ), is independent of the halo aggregation 
history (see Sec. [4} and so are also both the anisotropy and 
velocity dispersion profiles (SSMG) and the mass density 
profiles (SVMS) setting the truncation radii, the truncated 
subhalo number density profiles, N tT (r, Mg r ), must be also 
independent of the halo aggregation history. In other words, 
all the properties derived so far should not depend on the 
halo aggregation history. This would explain why all sub- 
structure properties derived from numerical simulations re- 
porting to haloes with very different aggregation histories, 
show such small scatters. 



6 THE EFFECTS OF DYNAMICAL FRICTION 

But things are not that simple. The spatial distribution 
of subhaloes is also affected by dynamical friction result- 
ing from gravitational two-body interactions between sub- 
haloes themselves and between subhaloes and diffuse dark 
matter particles. As in a major merger, the radial location 
of subhaloes suffers an important scrambling, the effects of 
dynamical friction that have previously taken place are es- 
sentially erased. Therefore, the effects of dynamical friction 
depend on the time elapsed since the last major merger. 
This means that, contrarily to all the processes previously 
mentioned, dynamical friction may lead to significant differ- 
ences between haloes according to their aggregation history 
(through the time elapsed since the last major merger). 

The large number of tiny subhaloes and of diffuse dark 
matter (in numerical simulations) suggests that dynamical 
friction should be very effective, at least in the case of the 
most massive subhaloes more prone to suffer it. But, this 
should have repercussions on the smooth structure and kine- 
matics of haloes, which should then depend on the halo ag- 
gregation history, while there is no clear sign of such a de- 
pendence in simulations (see SVMS and references therein). 
Furthermore, as a consequence of dynamical friction, the 
most massive subhaloes should lie closer to the halo cen- 
tre than less massive ones, whereas subhaloes of different 
masses show identical scaled number density profiles. And 
the only minor difference is rather of the opposite sign: the 
more massive the subhalo, the larger the mi nimum radius 
reach ed by its scaled number density profile (|Angulo et al.l 
2009). In particular, there is no sign of very massive haloes 
being accumulated at the halo centre. Therefore, simulated 
haloes show no apparent effect of dynamical friction. 

The only way to escape this paradox seems to be that 
the effects of dynamical friction may be present but go un- 
noticed. Does this make sense? If the only subhaloes having 
had time to suffer significant dynamical friction since the last 
major merger were the most massive ones, it could be very 
difficult to detect it because the number of those subhaloes 
is so small that their number density profile is quite uncer- 
tain (due to large Poisson errors). Certainly, they could be 
quite numerous at the halo centre where they should tend to 

7 The radial mapping of haloes is only preserved during accre- 
tion periods. However, the structural and kinematic properties of 
haloes resulting from a major merger arc indistinguishable from 
those shown by haloes grown by PA (see SVMS and SSMG). 



© 0000 RAS, MNRAS 000, 000-000 



12 Salvador-Sole, Serra, & Manrique 



accumulate, but, when a very massive subhalo falls to the 
halo centre and merges with the massive subhalo already 
lying there or simply settles down well-centred, it becomes 
invisible because it mimics the central part of the host halo. 
This effect could still manifest itself through an increased 
amount of small subhaloes near the halo centre, correspond- 
ing to old second-level subhaloes (in the massive subhaloes 
having fallen to the halo centre) converted to first-level ones. 
But such an effect should only be observable in case (a), be- 
cause, in case (b), subhaloes are mainly made of diffuse dark 
matter. In this sense, the lack of such an indirect proof of dy- 
namical friction would be an additional argument in favour 
of case (b) when trying to model the results of A-body sim- 
ulations. 



7 SUMMARY AND CONCLUSIONS 

The evolution of nested peaks in the filtering of the pri- 
mordial density field truthfully traces the evolution of halo 
substructure developing at all levels as a consequence of halo 
accretion and major mergers. Thus, the peak formalism can 
be used to describe subhalo abundance in typical haloes. 
Moreover, taking into account that haloes growing by PA 
develop from the inside out, it can be used to derive the 
non-truncated subhalo number density profiles per infinites- 
imal mass for subhaloes of different masses and, making use 
of the halo structural and kinematic profiles modelled in 
SVMS and SSMG, one can correct those quantities from 
tidal truncation as well. 

The subhalo properties predicted in the ACDM cosmol- 
ogy for Milky Way mass haloes are in very good agreement 
with those found in numerical simulations provided dark 
matter within subhaloes is essentially in the form of diffuse 
particles. The only slight deviations found, in this case, be- 
tween the theoretical predictions and the results of numer- 
ical simulations seem to be due to the non-fully accurate 
(sub)halo M-c relation used. More accurate M-c relations 
drawn e.g. from the SVMS model of halo structure would 
be welcome. 

But the true subhalo properties expected on pure the- 
oretical grounds in a real CDM universe rather correspond 
to those predicted under the opposite assumption that all 
dark matter in subhaloes is locked in higher-order level sub- 
haloes. Accurate predictions are also given for this more 
realistic case. The most striking result is that there should 
be, in this case, two orders of magnitude more subhaloes 
than usually thought on the basis of the results of A-body 
simulations. This might have important implications for the 
detectability of CDM, but it does not affect the dwarf galaxy 
abundance estimated from CDM simulations. 

In any of these two scenarios, subhalo properties are 
expected to be independent of the halo aggregation history. 
This means that, despite having been derived under the PA 
condition, all the previous quantities should hold for haloes 
having suffered major mergers. The only effects that could 
depend on the halo aggregation history are those due to dy- 
namical friction. According to the present results and those 
found in SVMS and SSMG, dynamical friction can be ne- 
glected as long as we are interested in modelling haloes in 
current simulations. However, it might have visible effects 
in the real CDM universe as well as in future, higher res- 



olution, A-body simulations. Of course, dynamical friction 
is also expected to have important consequences in baryon 
physics ignored in the present study. 
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APPENDIX A: TRUNCATED SUBHALO NUMBER DENSITY IN CASE (a) 

In case (a), i.e. M™ 111 = Ms, the subhalo number density per infinitesimal mass corrected for truncation (eq. [28]) takes the 
form 

f M _ _ fR(Ms) 

N tr (r,M s ) = N(r,M s ) + dM s N(r,M s ) _ dfA^ (f,M s ). (Al) 

By partial integration, equation (|A1|) leads to 

N(r,M s ) = N(r,M s )+ dM s N(r, > M s )-^- / df iV£ (f, M s ) . (A2) 

Jm s dM s y fltr(r ,M s) Ms ' M(r) 

Taking into account that both expressions (1A1[) and (|A2f) hold for any arbitrary value of Ms, we are led tcjfl 

a l ^ l r H ( M s) 

» [Ar(r ,<Ms)] = -^ / _ dFJV| (f,Jf B ), (A3) 

dMs dMs JR^(r,M s ) b W 

implying 

pR(M s ) 

3l(r,M s )N(r,<M s )= _ dftfg (f.Ms), (A4) 

where the function K(r, Ms) is the unknown integration constant. Choosing Ms equal to M(r) and taking into account that 
haloes grow inside-out, the double subindex "Ms,M(r)" in the truncated subhalo number density in the integrand on the 
right can be chosen equal to "Mo, Mo" without any loss of generality (see the meaning of such a double subindex in eq. |28|). 
That is, for that particular value of Ms, the subhalo is a clone of the host halo, except for the fact that it has not grown 
since it was aggregated by the host halo. In particular, its (sub)subhalo number density corrected for truncation is identical 
to that of the host halo itself (the subindex "Mo, Mo" can be omitted) and the (sub)subhalo with mass Ms is found at the 
same minimum radius r(Ms) as in the host halo. Consequently, equation (|A4[I takes the form 



U(r,M s )N(r)= dfN tz (f,M s ), (A5) 

Jmaijr(M s ),S' r [r,Jf(r)]} 

where we have taken into account that the N[r, < M(r)] is but the total subhalo number density at r, denoted as N(r). 
At small r, we have r(Ms) ^ i? tr [r, M(r)] and differentiation of equation (|A5[) leads to 

N tr (r,M s ) = [n(r,M s )N(r)} . (A6) 
dr 

Substituting N (r, Ms) given by equation (|A6f) into equation (|A1|I . taking into account equation (|A4I) and the partial inte- 
gration of the Ms-integral in the resulting expression, we arrive at 



[K(r, M s )N(r)] = N(r, Ms) + ~K(r, M s )N{r) N{r, > M s ) 



2 _N(r,>M s ) 



N(r) 



(A7) 



For any reasonable (large enough) value of Ms, N(r, > Ms)/N(r) is negligible in front of unity, so equation (|A7|) takes the 
simple form 

[3?(r, M s )N(r)] = N(r, M s ) + [5ft(r, M s )N(r)] N(r, > M s ) . (A8) 

dr 

This is a differential equation for K(r, Ms), which can be solved for the initial condition 5R[r(Ms), Ms] = implied by equation 
(|A5|l , Then, replacing the solution K(r, Ms) into equation (|A6[) . we obtain the wanted number density per infinitesimal mass 
of truncated subhaloes, N tT (r,Ms). 

At a large enough r, hereafter denoted by r cn d, the condition r(Ms) < R tr [r, M(r)] will be finally met and this solution 
will no longer hold. In this new regime, differentiation of equation (|A5|) leads to 



AT tr (r,M s ) = j- [$l(r,Ms)N(r)] - ^I^MMl N tr {R u '[r, M(r)], M s } . (A9) 



8 This is a physical rather than mathematical implication. The dependence on Ms in the integrands does not allow one to strictly prove 
the equality. But only a very unlikely conspiracy would make it possible to balance any arbitrary change in the integration limits by that 
produced in the integrands if they were not really equal. 
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Substituting N (r, Ms) given by equation (1A9() into equation (IA1|) . taking into account equation ()A4[1 and integrating by 
parts the integral over Ms in the resulting expression, we obtain 



dr 



(r,M s )N(r) 



dR u [r, M(R)] 
dr 



7V tr {7? tr [r,M(r)],Ms} 



= N(r, M s ) + - SR(r, Ms) iV(r) JV(r, > Ms) 



iV(r, > M s ) 
N(r) 



(A10) 



which, for any reasonable (large enough) value of Ms, reduces to 

A [M(r,M s )N(r)] + ^IhMi^l N *{ R *[ r , M(r)], Af s } 

= AT(r, M s ) + [5R(r, Ms)JV(r)] iV(r, > M s ) . (All) 

As ii tr [r, M(r)] is smaller than r, the function AT tr {ii tr [r, M(r)], Ms} has been previously obtained in the range of small 
r, so the differential equation (|A11|) can then also be solved for the function 5R(r, Ms) with the initial condition given by 
the value of 5R(r, Ms) at r cn( j. Once 5?(r, Ms) has been determined, we can replace it in equation (|A9|) to obtain the wanted 
function N tT (r, Ms) in the new radial range. 

In fact, given that r(Ms) is greater than R tr {R tT (r, Ms), M[R tr (r, Ms)]} in the relevant subhalo mass range (i.e. except 
for Ms < 10 3 Mq), the differential equation (|A11[) can be solved analytically. Indeed, equation (|A4|) for r = R tT (r, Ms) then 
takes the form 



»[ J R tr (r,Ms),M s )]A(r) = 



R tr {r,M s ) 



drN tT (f,M s ). 



(A12) 



r(M s ) 



Thus, by differentiating it with respect to R tT (r, Ms) and replacing the resulting expression for N tr {R tT [r, M(r)], Afs} into 
equation (|A11|) . we arrive at 



din 
dlnr 



(r,M s )N(r)}=ri 1 



■N(r,Ms)/[fH(r,M 3 )N(r)] 
1 + N(r, > M s ) 



(A13) 



N(r, Ms) is much smaller than N(r) and J?f(r, Ms) is much greater than one in this radial range, except for Ms > 10 10 Mq, as 
is checked a posteriori from equation (|A5j) . Consequently, for Ms < 10 10 M© we can neglect the term N(r, Ms)/[3?(r, Ms)N(r)] 
on the right of equation (|A13[l , which leads to the following quite accurate solution 



R(r,M s )JV(r) = i4(fw) exp 



df 



N(f, > M s ) 
1 + N(f,> Ms) 



(A14) 



with A(r cn d) an integration constant whose value is obtained by continuity of the solution at r cll d. Finally, differentiating 
equation (|A14[) and taking into account equation (|A4[) we are led to 



xrtr/ v ^(rend) AT(r, > Mg) 

^ (r ' Ms)= l + iV(r,>M s ) 6XP 



df 



N(r, > M s ) 



l + AT(f,> M s ) 



(A15) 
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